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ABSTRACT 

. We have studied the submillimetre (submm) properties of the following classes of near- 

er ' infrared (NIR)-selected massive galaxies at high redshifts: BzisT-selected star-forming 

galaxies (BzKs); distant red galaxies (DRGs); and extremely red objects (EROs). We 
used the SCUBA HAlf Degree Extragalactic Survey (SHADES), the largest uniform 
submm survey to date. Partial overlap of SIRIUS/NIR images and SHADES in SXDF 
has allowed us to identify 4 submm-bright NIR-selected galaxies, which are detected in 
the mid-infrared, 24 fim, and the radio, 1.4 GHz. We find that all of our submm-bright 
NIR-selected galaxies satisfy the BzK selection criteria, i.e. BzK = (z — K)ab — (B — 
z)ab > —0.2, except for one galaxy whose B — z and z — K colours are however close 
to the BzK colour boundary. Two of the submm-bright NIR-selected galaxies satisfy 
all of the selection criteria we considered, i.e. they belong to the Bzif-DRG-ERO 
overlapping population, or 'extremely red' BzKs. Although these extremely red BzKs 
are rare (0.25 arcmin -2 ), up to 20% of this population could be submm galaxies. This 
fraction is significantly higher than that found for other galaxy populations studied 
here. Via a stacking analysis, we have detected the 850-/im flux of submm-faint BzKs 
and EROs in our SCUBA maps. While the contribution of z ~ 2 BzKs to the submm 
background is about 10-15% and similar to that from EROs typically at z ~ 1, BzKs 
have a higher fraction (~ 30 %) of submm flux in resolved sources compared with 
EROs and submm sources as a whole. From the SED fitting analysis for both submm- 
bright and submm-faint BzKs, we found no clear signature that submm-bright BzKs 
are experiencing a specifically luminous evolutionary phase, compared with submm- 
faint BzKs. An alternative explanation might be that submm-bright BzKs are more 
massive than submm-faint ones. 
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1 INTRODUCTION 

The first 3-4 billion years in the history of the universe ap- 
pears to have been a very important epoch for the forma- 
tion of present-day massive galaxies, in particular for ellipti- 
cal galaxies and bulges with masses of M > 10 11 Mq. Recent 
wide field near-infrared (NIR) surveys in conjunction with 
multi-band optical images or spectroscopy have resulted in 
a profusion of samples of massive galaxies at high redshifts 
(e.g. Franx et al. 2003; Fontana et al. 2003; Dickinson et al. 
2003; Fontana et al. 2004; Cimatti et al. 2004; Glazebrook 
et al. 2004; Daddi et al. 2004; Juneau et al. 2005; Labbe 
et al. 2005; Drory et al. 2005; Reddy et al. 2005; Kong et al. 
2006; Papovich et al. 2006; Daddi et al. 2007; Arnouts et al. 
2007; Pozzetti et al. 2007). At z ~ 2, the physical prop- 
erties of massive galaxies are very different from those of 
present-day ones, i.e. a large fraction of them are intensively 
star-forming galaxies (e.g. Daddi et al. 2005), rather than 
passively evolving, like nearby ones. A notable character- 
istic of these massive star-forming galaxies is that, in some 
cases, a large total stellar mass ~ 10 1 Mq could be formed in 
the time-scale of a single starburst event, ~ 0.1 Gyr (Perez- 
Gonzalez et al. 2005; Caputi et al. 2006; Papovich et al. 
2006). 

Massive, intensively star-forming galaxies have also 
been found in submillimetre (submm) surveys (Smail et al. 
1997; Hughes et al. 1998; Eales et al. 1999; Scott et al. 2002; 
Borys et al. 2003; Mortier et al. 2005; Coppin et al. 2006). 
These submm galaxies turn out to be typically starbursts 
or starburst/AGN composite systems at z ~ 2 (Chapman 
et al. 2005; Alexander et al. 2005, and references therein), 
which are massive (Swinbank et al. 2004; Greve et al. 2005). 
They are likely to become passively evolving massive galax- 
ies at lower redshifts (e.g. Smail et al. 2004; Blain et al. 
2004; Takagi et al. 2004). The star formation rates (SFRs) 
of submm-bright galaxies are estimated to be extraordinar- 
ily high, ~ 10 3 M Q yr _1 (e.g. Greve et al. 2005; Alexander 
et al. 2005; Chapman et al. 2005; Pope et al. 2006; Kovacs 
et al. 2006), enough to produce a massive elliptical galaxy 
(> 3 L*) within ~ 1 Gyr. Determining the nature of these 
star-forming galaxies at high redshifts may be the key to 
understanding the process of massive galaxy formation. 

Efficient methods for selecting massive star-forming 
galaxies at high redshifts with one or two optical/NIR 
colours have played a significant role in studying the prop- 
erties of massive galaxy populations. These include BzK- 
selected galaxies (selected in B — z and z — K colours) , dis- 
tant red galaxies (selected with J — A' colour), and extremely 
red objects (selected with R — K or / — A colour). Hereafter 
we call these galaxy populations NIR-selected galaxies or 
NIRGs. 

Daddi et al. (2004) proposed an effective method for se- 
lecting massive star-forming galaxies at 1.4<z<2.5 along 
with the selection of passively evolving galaxies at a sim- 
ilar redshift range in a colour-colour diagram of B — z vs 

2 — K. It is expected to be reddening independent, i.e. ap- 
plicable to heavily obscured galaxies like ultraluminous in- 
frared galaxies (ULIRGs). Their luminosities at UV (red- 
dening corrected), infrared, X-ray and radio indicate that 
BzA-selected star-forming galaxies (BzKs) 1 with K < 20 

1 Since we are studying cross-identifications between BzK- 



are typically ULIRGs with SFRs of ~ 200 M yr _1 and red- 
dening of E(B - V) ~ 0.4 (Daddi et al. 2004, 2005; Reddy 
et al. 2005; Kong et al. 2006). 

In an alternative approach, Franx et al. (2003) em- 
ployed a colour cut of (J — A' s )v cg a > 2.3 to select distant 
red galaxies (DRGs). This colour criterion selects both pas- 
sively evolving and heavily obscured star-forming galaxies 
at 2 < z < 4.5 with strong Balmer or 4000A breaks. Recent 
Spitzer imaging at 24 ^m suggests that the bulk of DRGs 
are dusty star-forming galaxies (Webb et al. 2006; Papovich 
et al. 2006). Similar to BzA-selected star- forming galaxies, 
DRGs are likely to be ULIRGs at z ~ 2 (e.g. Knudsen et al. 
2005), but to also include a significant fraction of less lu- 
minous dusty star-forming galaxies at 1 < z < 2 (Grazian 
et al. 2006). Submm galaxies, another sample of ULIRGs 
at z ~ 2, could have a large overlap with BzKs or DRGs, 
but with systematically higher SFR than the average (e.g. 
Dannerbauer et al. 2006). 

Extremely red objects (EROs) defined with R — K or 
I — K colours have received particular attention as possi- 
ble optical-NIR counterparts of submm galaxies (e.g. Smail 
et al. 1999; Webb et al. 2003; Pope et al. 2005), since submm 
galaxies are expected to be heavily obscured by dust and 
therefore red. However, it turns out that only 10-30 % of 
submm galaxies, whose optical counterparts are identified 
with associated radio sources, are EROs (Ivison et al. 2002; 
Webb et al. 2003; Borys et al. 2004; Smail et al. 2004). The 
median redshift of EROs with A Vcg a < 19.2 is z = 1.1 ± 0.2 
(Cimatti et al. 2002), and therefore lower than that of re- 
solved submm-bright galaxies, BzKs and DRGs. Neverthe- 
less, there is certainly some overlap between their redshift 
distributions. 

It is important to understand the the physical relation- 
ships between these high-z massive star-forming galaxies 
which are selected in different ways (e.g. see Reddy et al. 
2005). The surface density of BzKs, DRGs and EROs is 
similar to each other, about 1-2 arcmin -2 with Av og a < 20 
(Daddi et al. 2004; Cimatti et al. 2002; van Dokkum 
et al. 2003). A similar surface density has been found for 
submm galaxies with Ssso^m > 1 mjy (e.g. Blain et al. 1999), 
although the surface density of bright (Ssso/jm > 8 mjy) 
submm galaxies in 'blank sky' surveys is at least an order of 
magnitude smaller. Understanding the physical origin of the 
overlap between submm galaxies and NIR-selected galaxies 
is not straightforward, given the diverse optical-NIR proper- 
ties of submm galaxies, which include those similar to less- 
attenuated UV-selected galaxies as well as EROs (e.g. Smail 
et al. 2002, 2004). 

Here we study the submm properties of BzKs, DRGs 
and EROs by using the SCUBA HAlf Degree Extragalac- 
tic Survey (SHADES - Mortier et al. 2005), providing a 
larger sample of submm galaxies for investigating the rela- 
tion between submm galaxies and NIRGs than the study 
of Reddy et al. (2005) for galaxies in the GOODS-North 
field. SHADES has mapped two regions in the Lockman 
Hole and in the Subaru/ XMM-Newton deep field (SXDF). 
In this study, we focus on a 93arcmin 2 sub-region of the 



selected galaxies and submm galaxies, star-forming galaxies se- 
lected with this method are simply referred to as BzK-selected 
galaxies or BzKs. 
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SXDF field of SHADES, in which both optical (BVRi'z') 
and N1R (JHK S ) photometry are already available. We de- 
scribe the data in Section 2 and the sample of NIR-selected 
galaxies in Section 3. In Section 4, we identify submm-bright 
NIR-selected galaxies and assess the overlap between submm 
galaxies and NIR-selected galaxies. In Section 5, we present 
the results of a statistical detection of NIR-selected galax- 
ies in the SCUBA map using a stacking analysis. In Sec- 
tion 6, we analyse the spectral energy distributions (SEDs) 
of submm-bright NIR-selected galaxies and also an average 
SED of submm-faint NIRGs, in order to investigate their 
physical properties. Summary of our study is given in Sec- 
tion 7. Throughout this paper, we adopt the cosmology of 
fi m = 0.3, Oa = 0.7 and H = 70kmsec _1 Mpc" 1 . All the 
magnitudes are given in the AB system, unless otherwise 
noted. 

2 THE DATA 

2.1 Optical/NIR data in the SXDF 

We have used the optical and NIR imaging data of the 
SXDF centred at [2 h 18 m 00 s , -5°00'00" (J2000)] published 
by Miyazaki et al. (2003). The NIR data were obtained with 
the University of Hawaii 2.2 m telescope with the Simultane- 
ous 3-colour InfraRed Imager for Unbiased Survey (SIRIUS; 
Nagayama et al. 2003). The final image covers an area of 
114arcmin 2 in the H- and A a -bands. For the J-band, we 
obtained data only for 77arcmin 2 , since a quarter of the 
field of view was not operational. The 5 a detection lim- 
its of the final NIR image are J = 22.8, H = 22.5 and 
A' s = 22.1 mag, through a 2" diameter aperture. From the 
final image in the A s -band, we have detected 1308 objects 
having K s < 22.1 mag using SExtractor (Bertin & Arnouts 
1996). Throughout this work we use the total magnitudes 
in the optical-NIR bands measured with the MAG_AUTO 
algorithm in SExtractor. 

For these A s -band detected objects, we made a multi- 
band catalogue using optical images of Subaru/Suprime-cam 
(Miyazaki et al. 2002) in B, V, R, i', z bands 2 . The limiting 
magnitudes (5 ex) of our images are B — 27.1, V = 25.9, 
R = 26.3, i' — 25.7, z = 25.0 mag through a 2" diameter 
aperture. The FWHM of the point-spread function is 0.98" 
in both the optical and NIR images. Stars were identified 
with the colour criterion of B — K B < 1.583(£> — i') — 0.5 and 
FWHM < 1.2" (Miyazaki et al. 2003), and removed from 
the analysis. 

2.2 Submm/radio data from SHADES 

Details of the survey design and observing strategy with 
JCMT/SCUBA are given in Mortier et al. (2005). Here we 
describe the SHADES survey only briefly. Using the SCUBA 
instrument on the JCMT we have obtained jiggle maps of 
the SXDF and also the Lockman Hole in grade 2-3 weather 
(rcso = 0.05 - 0.1). The FWHM of the SCUBA beam 
is 14.7" and 7.8" at 850 and 450 /mi, respectively. The 
SHADES observations with SCUBA were continued from 

2 Our data were obtained during the commissioning phase of the 
Suprime-cam. 



December 2002 to the decommissioning of SCUBA in June 
2005. Here we use the SCUBA data acquired until 1st Feb, 
2004. 

The SCUBA data have been reduced by four indepen- 
dent groups in the SHADES consortium. Practical methods 
of data reduction, such as flux calibration, extinction cor- 
rection, and map making depends on the reduction group 
and are summarised in Coppin et al. (2006). The extracted 
sources from four different SCUBA maps have been com- 
pared and then combined to produce the SHADES source 
catalogue, which is expected to be the most reliable source 
catalogue from SCUBA surveys. In our analyses, we focus 
on two particular SCUBA maps out of four, along with the 
SHADES source catalogue, in order to make the sample of 
submm-bright NIRGs as complete as possible. 

We have chosen to use only regions with noise less 
than 3 mjy at 850 /im. The total area of this region is 
~250 arcmin 2 . The field centre for the SHADES map of the 
SXDF, 02 h 17 m 57 s .5, -05°00'18".5 (J2000), is offset from 
that of the SIRIUS observations and hence our sub-mm 
maps do not cover the whole region covered by the NIR data. 
In one of the four SCUBA maps ('Reduction B' in Coppin 
et al.), this overlap is 93 arcmin 2 , with median noise levels 
of 2.0 mjy and 18.4 mjy at 850 and 450 /im, respectively. 

Wide- field 1.4-GHz radio images of SXDF were ob- 
tained using the VLA during 2003. Around 60 hr of inte- 
gration were salvaged, following a prolonged failure of the 
correlator, comprising data from the A, B and C configu- 
rations, with an approximate 9:3:1 ratio of recorded visi- 
bilities, evenly distributed in three pointings separated by 
15 arcmin. The final images were mosaiced together, after 
correcting for the response of the primary beam. The re- 
sulting noise level is around 7 /ijy beam -1 in the best regions 
of the map (rather higher near bright, complex radio emit- 
ters) with a synthesised beam measuring 1.87 x 1.65arcsec 
(FWHM), major axis 22° east of north. The data and their 
reduction are described in detail in Ivison et al. (2007). 



2.3 Spitzer data from SWIRE 

The whole area of our A s -band image is covered by the 
Spitzer Wide-area InfraRed Extragalactic survey (SWIRE; 
Lonsdale et al. 2003). We used the IRAC+24/mi cata- 
logue from version 2.0 data products (released in Sum- 
mer 2005), available from the NASA/IPAC Infrared Ser- 
vice Archive (IRSA). We adopted the Kron fluxes, which 
are MAG_AUTO fluxes from SExtractor, in order to com- 
pare with the ground-based optical-NIR photometry. This 
catalogue includes IRAC sources which are detected at both 
3.6 /im (> 10 <t) and 4.5 /im (> 5a). The astrometric error 
of the IRAC sources is small enough to easily identify coun- 
terparts of NIRGs with the angular separation of < 0.5". In 
the SIRIUS/ A s -band region, the faintest sources at 24 /im 
in the catalogue have a flux of ~300/iJy. Cross-correlation 
between the SHADES sources and the SWIRE sources in 
SXDF is given in Clements et al. (2007, in preparation). 
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3 SAMPLE OF NIR-SELECTED GALAXIES AT 
HIGH REDSHIFTS 

3.1 EROs and DRGs 

Following Miyazaki et al. (2003), we define EROs as objects 
with R — K B > 3.35, which is equivalent to (R — K )v cg a > 5, 
i.e. the widely used criterion of EROs. In order to obtain 
the DRG sample, we use the threshold of J — K s > 1.32 (i.e. 
>2.3 in Vega magnitudes). Within the SHADES area, we 
have detected 201 EROs and 67 DRGs in total, among which 
39 objects satisfy both R — K s > 3.35 and J — K s > 1.32. 
In Figure 1, we show a colour-colour plot with R — K B and 
J — K B for its-band detected objects, along with the adopted 
colour criteria. The statistics of NIRGs are summarised in 
Table 1. 

We derived a surface density of 2.18 ± 0.14 and 1.09 ± 
0.12 arcmin -2 (Poissonian errors) for EROs and DRGs, re- 
spectively. By using the same data, Miyazaki et al. (2003) 
found a good agreement in the surface density of EROs 
with those of the other surveys with areas of > 50 arcmin 2 
(Cimatti et al. 2002; Daddi et al. 2000; Thompson et al. 
1999). The surface density of DRGs is coincidentally very 
similar to that derived by van Dokkum et al. (2003), al- 
though the A" s -band limiting magnitude in van Dokkum 
et al. (2003) is ~ 0.7 mag deeper than ours. Thus, we find a 
higher surface density of K B < 20 DRGs than in van Dokkum 
et al. (2003). This could be partly because of the contam- 
ination of the sample at K B > 21.5 where the detection of 
DRGs at J-band is possible only below 4a. 



3.2 BzK-selected star-forming galaxies 

Daddi et al. (2004) proposed a joint selection of star-forming 
galaxies and passively evolving galaxies in the redshift range 
of z = 1.4-2.5 in the z — K vs B — z diagram (here- 
after l BzK diagram'). Following Daddi et al. (2004), we 
choose BzA'-selected galaxies with BzK > —0.2, where 
BzK = {z — K) — {B — z). In order to adjust the selec- 
tion in our photometric bands to that of Daddi et al., we 
compared the stellar sequence in the BzK diagram with 
that of Daddi et al. (2004). To match the stellar sequence, 
we applied the following colour corrections: (B — z)Daddi = 
(B - z) + 0.2 and (z - K) D addi=(z' - K) - 0.2. After this 
correction, we obtained 132 BzA-selected galaxies within 
the SHADES area. 

We obtained a surface density of 1.1 arcmin -2 for K s < 
20 (Vega). This surface density is consistent with that 
of Daddi et al. (2005), who found 169 BzKs within 154 
arcmin , including X-ray detected objects, in the GOODS 
North region. Thus, we assume there are no systematic dif- 
ferences between our sample and that of Daddi et al. (2005) . 
In Figure 2, we plot the raw BzK colour (i.e. with no colour 
corrections) of AVband detected objects against K B magni- 
tudes along with the photometric errors. 



4 IDENTIFICATION OF SUBMM-BRIGHT 
NIRGs 

4.1 Method of identification 

In order not to miss any candidate submm-bright NIRGs, 
we utilized two SCUBA 850-^tm maps which have different 
pixel size produced by independent reductions and reported 
in Coppin et al. (2006 - Reductions B and D). The pixel sizes 
of the SCUBA 850-^mi map from Reductions B and D are 
1" and 3", respectively. One benefit of using SCUBA maps, 
rather than the SHADES source catalogue, is the possibility 
of identifying additional submm sources which are just be- 
low the limit of the SHADES catalogue. This is a reasonable 
approach, since we know the sky positions of targets before- 
hand (and see Appendix A for a statistical discussion). 

We adopted the threshold of S/N > 3 for the detection 
of submm fluxes at the positions of NIRGs. In order to con- 
firm the detection of NIRGs in SCUBA maps, we need to ex- 
clude the possibility of; 1) chance association with a nearby 
SCUBA source; and 2) false positive detections. Sources in 
the SHADES catalogue should be quite reliable, since they 
are extracted using 4 independent reductions. If there are 
no SHADES catalogue sources corresponding to detected 
NIRGs, we need to pay special attention to check the reli- 
ability of the detection in the SCUBA map. To overcome 
these problems, we use the Spitzer images at 24 and 
the VLA radio images at 1.4 GHz as they have proved to 
be useful in previous studies (e.g. Ivison et al. 2002; Egami 
et al. 2004) . The identification of SHADES sources has been 
undertaken by Ivison et al. (2007). We confirmed that our 
results are consistent with their identifications. For non- 
SHADES sources, we require detection at 24 fim or in the 
radio for secure identification of submm-bright NIRGs. 

We also obtained 450-/^m images simultaneously with 
the 850-^mi images, although the weather conditions allo- 
cated for SHADES were not good enough to detect typical 
850-/im sources at 450 /im in the same integration time. We 
have however used the 450 fim map to constrain the 450- 
fim flux for the detected 850-/^m sources. The reduction of 
450 [im data is fully described in Coppin et al. (2006). 



4.2 Results 

Among 307 NIRGs in the SHADES/SIRIUS AVband area 
(245 fall in J-band area), we detected 5 NIRGs at 850/im 
as a result of photometry at the position of NIRGs in the 
SHADES B-map (i.e. the SCUBA 850 /im map from Re- 
duction B). We also detected 5 NIRGs at 850/j,m in the 
SHADES D-map, four of which are common with those de- 
tected in the SHADES B-map. The results of the sub-mm 
detections from both SHADES maps are summarised in Ta- 
ble 2. All of the submm-detected NIRGs except fo ID1390 
have corresponding SHADES catalogue sources. Although a 
submm source for ID1390 (SXDF850.62) was included in a 
preliminary SHADES source catalogue, this source did not 
satisfy the final significance limit adopted for the SHADES 
catalogue presented by Coppin et al. (2006). However, we 
found that ID 1390 is detected both at 24 fim and in the ra- 
dio, and therefore we conclude that this is a real submm 
galaxy (see also Appendix A). In the 450-/im map, we find 
no convincing detections at the positions of NIRGs. This 



Submm properties of NIR-selected galaxies in SHADES/SXDF 5 



is not unexpected given the high noise level in the 450/ira 
maps. 

Three (ID300, 445, 912) out of six submm-detected 
NIRGs are identified as optical counterparts of correspond- 
ing SHADES sources in Ivison et al. (2007). Additionally, we 
find ID1390 as an optical-NIR counterpart of SXDF850.62 (a 
non-SHADES source), since this galaxy is detected at 24 fim 
and in the radio as noted above. In summary, we identify 
four submm-bright NIRGs. Figure 3 shows the 850-/im con- 
tours plotted over K B -h&n& images of submm-bright NIRGs. 
Figure 3 also shows the radio contours over i?-band images 
of the same NIRGs. The remaining two submm-detected 
NIRGs (ID718 and ID1133) are detected neither in the radio 
nor at 24 /jm. Such chance associations are expected and the 
number is consistent with the estimate given in Appendix 
A. For a submm source associated with ID1133, Ivison et 
al. (2007) found another reliable optical counterpart with a 
radio detection. Considering these results, we exclude both 
ID718 and 1133 from further analyses. 

In order to check additional candidates for optical-NIR 
counterparts of the SHADES catalogue sources, we also 
searched for NIR-selected galaxies within 7" radius of the 
SHADES positions. Although we found two additional as- 
sociations (ID475 and 579), they are not detected either at 
24 fim or in the radio. Therefore, this alternative method 
does not produce any further submm-bright NIRGs for our 
sample. 

4.3 Relation between submm galaxies and NIRGs 

In Table 3, we tabulate the optical-NIR properties of 
submm-bright NIRGs. Note that all of the submm-bright 
NIRGs are BzTf -selected galaxies, except for ID1390 with 
BzK = —0.36, although it satisfies BzK > —0.2 before the 
colour correction on B — z and z — K. This also means that 
no EROs and DRGs which are clearly non-BzK galaxies 
are detected at 850 /im. The distribution of submm-bright 
NIRGs in the B — z vs z — K colour-colour diagram is dis- 
cussed in detail in Section 6.1. 

Two galaxies, ID300 and ID445, satisfy all the colour 
selection criteria we adopted, i.e. BzK > —0.2, R — K > 
3.35 and J — K > 1.32. Hereafter, we refer to this BzK- 
DRG-ERO overlapping population as extremely red BzK- 
selected galaxies. We found only 16 extremely red BzKs in 
the combined K s and SHADES areas, corresponding to a 
surface density of only 0.25 ± 0.05 arcmin -2 . This means 
that 12 ± 8% (2/16) 3 of extremely red BzKs are submm 
galaxies. Although the sample is very small, nevertheless 
it may indicate that this rare NIR galaxy population has 
a much higher fraction of submm galaxies, compared with 
the other optical-NIR selected galaxy populations already 
studied. For example, the fractions of submm galaxies which 
are BzKs, EROs and DRGs in our sample are 4/132 (3 ± 2 
%), 3/201 (1.5 ± 1 %) and 3/67 (5 ± 3 %), respectively. 

We found that there are 20 (13) SHADES sources in 
the SIRIUS/ifs-band (J-band) area, among which three are 



3 The errors on fractions, /, are given by y ^ where n is 

the total number of the sample. This gives a rough estimate based 
on the normal approximation to the binomial distribution. 



BzK-selected and two satisfy both the ERO and DRG se- 
lection. The resulting fractions of BzKs, EROs and DRGs 
in our sub-sample of SHADES sources are 15 ± 8 %, 10 ± 7 
% and 15 ± 10 %, respectively. Thirteen SHADES sources 
out of 20 have a robust radio identification (Ivison et al. 
2007). If we confine our attention to SHADES sources with 
a robust radio identification, the above fractions of NIRGs 
are higher by a factor of ~ 1.5. 

Ivison et al. (2002) found an ERO fraction of 33 % 
(6/18) in SCUBA sources with robust radio identification 
from the 8-mJy survey, in which one ERO is too faint to 
be detected at our detection limit in the _K" s -band. Reddy 
et al. (2005) studied the BzK and DRG fraction of submm 
galaxies with Ssso/jm ^ 5 mjy in the GOODS-North field. 
Out of 11 radio-detected submm galaxies, they found that 
5 (45±15%) and 3 (27 ± 13%) objects satisfy the BzK and 
DRG criteria, respectively. For these submm-bright BzK 
and DRG samples, 2 objects each are too faint to be de- 
tected at our detection limits. Subtracting these Jf-band 
faint sources would give the BzK and DRG fractions of 
27±13 % and 9±8 %, respectively. On the other hand, Daddi 
et al. (2005) found only one _Bz7f-selected galaxy in the list 
of submm sources in Pope et al. (2005) in the GOODS-North 
field. These results are not inconsistent with ours, consider- 
ing the large errors, due to a small sample size and cosmic 
variance. 

Since the redshift distribution of BzKs has a large over- 
lap with that of submm galaxies, and both classes of object 
are star-forming galaxies, one could expect a large fraction 
of BzKs within the submm galaxy population. We roughly 
estimate the expected fraction of BzKs in submm galaxy 
samples as follows, based on the results of the redshift survey 
of radio-detected submm galaxies by Chapman et al. (2005) . 
In their sample of radio-detected submm galaxies, the frac- 
tion of galaxies at 1.4 < z < 2.5 is about 45 %. This fraction 
should be corrected for the incompleteness of the sample, 
mainly due to the redshift desert at 1.2 < z < 1.8. Given 
the 20% incompleteness at z ~ 1.2-1.8 estimated by Chap- 
man et al. (2005), the fraction of galaxies at 1.4 < z < 2.5 
may be in the range of 35-45 %. If BzKs are an almost 
complete sample of star-forming galaxies at 1.4 < z < 2.5, 
including heavily obscured objects, most of radio-detected 
submm galaxies at 1.4 < z < 2.5 may be selected as BzKs. 
In our SHADES/SIRIUS field, there are 13 radio-detected 
SHADES catalogue sources. Therefore, one may expect that 
4-6 of them are BzKs. From the K-band magnitude of 
radio-detected submm galaxies of Smail et al. (2004), we 
found that ~65 % of radio-detected submm galaxies are de- 
tectable at the detection limit of our _R" s -band image. Thus, 
we expect that at most ~4 radio-detected submm galaxies 
would be BzKs in our sample. This rough estimate is con- 
sistent with the results above, since we found 3 such radio- 
detected submm galaxies. Although the sample is small, this 
suggests that the BzK selection technique is effective even 
for submm galaxies, i.e. heavily obscured galaxies with ex- 
tremely high SFR. Thus, submm galaxies at z ~ 2 may 
have the largest overlap with BzK-selected galaxies, when 
compared to EROs and DRGs. In Section 6.2.3, we discuss 
the implications of this, using a radiative transfer model of 
starburst galaxies. 

The results presented here are hampered by small sam- 
ple size, partly due to the limited overlap between our NIR 
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images and the SHADES field. Both of the SHADES fields 
have been covered by the UKIRT Infrared Deep Sky Survey 
(UKIDSS - Lawrence et al. 2006). In the near future it will 
therefore be possible to study the submm properties of a 
large sample of NIRGs using the UKIDSS data. 



5 STACKING ANALYSIS FOR SUBMM- FAINT 
NIRGs 

We can try to extract a little more information about the 
submm galaxies/NIRG overlap by looking at the statistical 
correlations within the images. Here we estimate the average 
850- /im flux of submm-faint NIRGs and the contribution to 
the extragalactic background light (EBL) from each class of 
objects. 

We first eliminated the effects of resolved submm 
sources in our SCUBA maps. In the B-map, we excluded 
regions within 7" of the SHADES catalogue sources. In the 
D-map, all the SHADES sources and their associated nega- 
tive off-beams were removed before measuring the flux. For 
this cleaned D-map, we confirmed that random positions 
give a variance-weighted average flux of zero. Although a 
non-SHADES source for ID 1390 is not removed from the 
maps, its effect is negligible. 

In Figures 4-6, we show the histogram of measured flux 
at the position of NIRGs in the B-map, compared with that 
of the map as a whole. The significance of the difference be- 
tween the distributions of the map and the measured fluxes 
was estimated with the Kolmogorov-Smirnov test. The flux 
distributions of BzKs and EROs were found to be different 
from that of the map as a whole, with a significance of 2-3 a 
(see Table 4). For DRGs, we found no significant detection 
from the Kolmogorov-Smirnov test using 56 objects. As a 
comparison, Knudsen et al. (2005) detected an average flux 
of 0.74 ± 0.24 mjy for DRGs (excluding discrete sources) 
with a smaller sample of 30, but with a lower noise SCUBA 
map. A possibly large contamination in our DRG sample, 
owing to the shallow J-band data, might cause no significant 
detection. 

From the B-map, we obtained an average flux of 0.52 ± 
0.19 and 0.53 ± 0.16 mjy for BzKs and EROs, respectively. 
We measured the average flux of data between the first and 
the third quartiles of the sample, which is a robust measure 
against outliers. The errors are given by the median absolute 
deviation. Table 4 gives a summary of the stacking analysis 
using the B-map. The noise-weighted average flux of BzKs, 
EROs, and DRGs measured in the D-map is 0.64 ± 0.16, 
0.50 ±0.13 and 0.42 ±0.23 mjy, respectively, which are con- 
sistent with the results from the B-map. The derived aver- 
age fluxes of BzKs and EROs are consistent with previous 
studies in other sky regions (Webb et al. 2004; Daddi et al. 
2005) 4 . In the following discussion, we use the average flux 
of submm-faint BzKs and EROs from the B-map. 

We next estimate the contribution from each class of 
object to the EBL at 850 /im. The total flux from individ- 
ually detected sources is added to the estimated total flux 

4 Note that Daddi et al. (2005) derived an average flux of 1.0 ± 
0.2 mjy for '24/jm-detected' BzKs. Including 24 ^tm-undetected 
objects, this average becomes 0.63 ± 0.17mjy in the GOODS- 
North field 



from undetected objects. The values obtained for the EBL 
from BzKs and EROs are 3.8 ± 1.2 and 5.1 ± 1.5 Jydeg" 2 , 
respectively. The measured 850 pm EBL is 31Jydeg~ 2 in 
Puget et al. (1996) and 44Jydeg~ 2 in Fixsen et al. (1998), 
i.e. there is a relatively large uncertainty on the absolute 
level. The contribution from BzKs and EROs are both 10- 
15 % each, depending on the adopted value of the EBL. This 
modest contribution is also suggested by a comparison of 
the surface density of objects. For example, the surface den- 
sity of BzKs, ~ 1.5 arcmin~ 2 is a factor of ~ 5 smaller than 
that of >0.5 mjy SCUBA sources, ~ 7 arcmin -2 (Blain et al. 
1999; Cowie et al. 2002) at which flux density level the back- 
ground is close to complete. 

We also estimate the resolved fraction of the EBL orig- 
inating from a given galaxy population only. The resolved 
fraction of the EBL from BzKs is at least ~ 30 %, which 
could be higher than that found in EROs. The detected 
sources in our sample correspond to > 3.5 a sources when 
we measure the peak flux in the SCUBA map with the noise 
level of ~ 2 mjy. From the 8-mJy survey, which has a similar 
noise level, Scott et al. (2002) found that > 3.5 a sources ac- 
count for ~ 10 % of the EBL. Thus, the submm flux of BzKs 
is apparently biased high; i.e. a large fraction of the submm 
fluxes from BzKs are found in resolved bright sources. 

Recently, Wang et al. (2006) suggested that the major- 
ity of the EBL at 850 )im originates from submm-undetected 
galaxies at z<1.5. On the other hand, the majority of 
submm-detected galaxies, i.e. resolved sources, lie at z > 1.5 
(Chapman et al. 2005). This suggests that the resolved frac- 
tion of the EBL varies as a function of redshift. If we study 
the EBL only from galaxy populations at z > 1.5 like BzKs, 
the contribution from submm-undetected galaxies to the 
EBL would be significantly reduced. This may explain why 
the resolved fraction of the EBL from z ~ 2 BzKs is higher 
than that from EROs typically at z ~ 1. 



6 SED ANALYSIS OF SUBMM-BRIGHT 
NIRGs 

In this section, we investigate the physical properties of 
submm-bright NIRGs from the observed SEDs. Firstly, we 
examine the distribution of submm-bright NIRGs in the 
B — z vs z — K colour-colour plot, i.e. BzK diagram. Specif- 
ically, we compare the properties of submm-bright NIRGs 
with those of 24 ^im-detected NIRGs. Secondly, we apply 
an evolutionary SED model of starbursts to submm-bright 
NIRGs for a more detailed study of each object. Also, we 
study the properties of submm-faint NIRGs from an aver- 
age SED, which are compared with those of submm-bright 
NIRGs. 

6.1 Submm-bright NIRGs in the BzK diagram 

The fraction of BzKs in a galaxy population would de- 
pend on its redshift distribution; i.e. low-z galaxy popula- 
tion would have a small fraction of BzKs. Figure 7 shows 
the BzK diagram for all the A's-detected objects in the 
SXDF/SIRIUS field. While NIRGs have a wide range of 
B — z, i.e. < (B — z) < 6, submm-bright NIRGs favour the 
colour ranged <(B -~z) < 3, with BzK> - 0.2. This is 
not the case for general 24 /im-detected objects, which have 
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0.5 < (B — z) < 4 and a large fraction of non-BzA's. Un- 
like 24 ^in-detected NIRGs, submm-bright NIRGs disfavour 
xion-Bz K EROs and are better matched to the B — z colours 
of BzKs. This is consistent with the redshift distribution of 
radio-detected submm galaxies, having a median of z ~ 2 
(Chapman et al. 2005), which overlaps considerably with 
that of BzKs. On the other hand, 24 pm-detected objects 
include a large number of z < 1.5 galaxies (Rowan- Robinson 
et al. 2005), and have a small BzK fraction. 

Note that 24 ^in-detected NIRGs include fairly blue 
BzKs with (B — z) < 1, unlike submm-bright ones. The 
correlation between B — z and reddening for BzA'-selected 
galaxies is discussed by Daddi et al. (2004), suggesting that 
E(B-V) = 0.25(B- z + 0.1) for the Calzetti extinction law 
(Calzetti et al. 2000). According to this relation, (B — z) <1 
corresponds to E(B — V) <0.3. On the other hand, we de- 
rive E(B — V) ~ 0.5 for submm-bright BzKs on average 
from (B - z) = 2.1. Reddy et al. (2005) show that UV- 
selected 'BX/BM' galaxies occupy a region in the BzK di- 
agram similar to BzA-selected galaxies with (B — 
Since these UV-selected galaxies are less obscured by dust 
than submm galaxies (e.g. Smail et al. 2004), it is expected 
that most submm-bright NIRGs may avoid the colour region 
of (B -z)<l. 

The other possibilities for blue B — z colours include the 
presence of AGN. We performed a cross-correlation between 
the 24 /^m-detected BzK sample and X-ray sources in the 
XMM serendipitous source catalogue 5 . As a result, we found 
that the bluest four 24 /^m-detected BzKs with (B — z) < 1 
and (z — K) < 1 are associated with X-ray sources with an 
angular distance of 9 < 2", while the submm-bright sam- 
ple have no such associations. Therefore, bluer colours of 
24 ^im-detected BzKs appear to be better explained by a 
contribution from AGN. 

In summary, the colours of submm-bright NIRGs, 
1 < (B — z) < 3, indicate that they are obscured star-forming 
galaxies at z>1.4, with no obvious contribution from 
AGN to the observed optical-NIR fluxes. The nature of 
submm-bright NIRGs is further discussed below using multi- 
wavelength data and a theoretical SED model. 

6.2 Comparison with SED models 

6.2.1 SED fitting method 

We analyse the SEDs of submm-bright NIRGs using an evo- 
lutionary SED model of starbursts of Takagi et al. (2003) 
which has previously been applied to submm galaxies in Tak- 
agi et al. (2004). In this model, the equations of radiative 
transfer are solved for a spherical geometry with centrally 
concentrated stars and homogeneously distributed dust. We 
use the same model templates as those used in Takagi et al. 
(2004). We hereafter refer to this model as the StarBUrst 
Radiative Transfer (SBURT) model. Here we extend the 
wavelength range of the SBURT model to the radio by as- 
suming the observed correlation between far-IR and radio 
emission (Condon 1992), with a = —0.75 and q = 2.35, 



5 The XMM-Newton Serendipitous Source Catalogue, version 
1.1.0, XMM-Newton Survey Science Centre Consortium, XMM- 
SSC, Leicester, UK (2004) 



where a is the spectral index of radio emission (S v oc v a ) 
and q defines the luminosity ratio of far-IR to radio emission 
at 1.49 GHz (see Condon 1992). 

The fitting parameters of the SBURT model are the 
redshift, starburst age and compactness of a starburst re- 
gion (O). The evolutionary time-scale of the starbursts to is 
assumed to be 0.1 Gyr, which specifies both the gas infall 
rate and the star formation rate. The compactness of star- 
bursts is defined by r — 0(M, / 10 9 Mq)^ [kpc], where r and 
M* are the radius and stellar mass of the starburst region, 
respectively. The dust model is chosen from the MW, LMC 
or SMC models taken from Takagi et al. (2003). Following 
the results of Takagi et al. (2004), we adopt a top-heavy ini- 
tial mass function (IMF) with a power law index of x = 1.10 
(the Salpeter IMF has x = 1.35) for submm galaxies, which 
is necessary to reproduce the colour-magnitude relation of 
present-day elliptical galaxies. The lower and upper mass 
limits of the adopted IMF are 0.1 and 60 Mq, respectively. 
This particular choice of the IMF does not affect the follow- 
ing results, except for the derived stellar masses. 

The best-fitting SED model is searched for using a \ 2 ~ 
minimization technique from the prepared set of SED mod- 
els. We used all the available fluxes, except for 24 fim and 
radio. This is because, in the SBURT model; 1) the con- 
tribution from an AGN is not taken into account; 2) the 
24 /im flux depends on the details of the dust model; and 
3) the radio flux is separately calculated by using the em- 
pirical relation. The adopted 850-pm fluxes and errors are 
'deboosted values' (i.e. accounting for the effects of flux 
boosting on a low S/N threshold sample) taken from Coppin 
et al. (2006) or calculated with the same deboosting algo- 
rithm. We adopted a minimum flux error of 5 %, considering 
the systematic uncertainty of photometry from the different 
type of instruments covering a wide range of wavelength. 
The upper limits on fluxes are taken into account in the 
fitting process, i.e. models exceeding the 5 a upper limits 
are rejected. At rest frame UV wavelengths, the SEDs of 
heavily obscured starbursts like submm galaxies could de- 
pend on the inhomogeneity of interstellar medium, since the 
effects of photon leakage may dominate the resulting SED 
(e.g. Takagi et al. 2003). Also, the absorption by the inter- 
galactic medium is important in the rest frame below 1216 A, 
which depends on a particular line of sight to each galaxy. 
Considering these uncertainties, we quadratically added an 
additional 20% error for data at rest-frame UV wavelengths 
below 4000 A. 



6. 2. 2 Results of SED fitting 

In Figure 8, we show the best-fitting models for submm- 
bright NIRGs. The original and radio-extended SBURT 
models are depicted with solid and dotted lines, respectively. 
The fitting and derived model parameters are summarised 
in Table 5. 

For ID912, we found that the value of minimum \ 2 re- 
duced drastically when we excluded //-band data from the 
fitting analysis. This may indicate that the //-band photom- 
etry could have a large systematic error. The upper limit in 
the J-band suggests that this galaxy could be very red in the 
observed NIR colours. The noise level in the //-band may 
be too high to detect this galaxy. We checked the images of 
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ID912 and found that photometry of this galaxy could be 
affected by a nearby object. Hence, we adopt the best fitting 
model for ID912 without i/-band data. 

We find a significant underestimate of 850 /im flux for 
ID445 which has a very red z' — K B colour, although the 
resulting \ value indicates that the best-fitting model is 
statistically acceptable. Since the predicted radio flux is also 
well below the observed flux, the best-fitting model might be 
rejected with more accurate measurements in the submm. 
For this source, we may particularly need a more compli- 
cated multi-component model of starbursts, in which young 
heavily obscured molecular clouds are treated separately 
(e.g. Silva et al. 1998). 

The SED fitting suggests that submm-bright NIRGs 
could lie at the typical redshift range of submm galaxies 
(Chapman et al. 2005). In order to show the fitting error on 
photometric redshifts z p hot, we show contour plots of A\ 2 
projected onto the age-redshift plane in Figure 9 for models 
with the best dust model for each galaxy. We note that the 
derived z p hot and the far-lR-radio relation reproduce the 
observed radio flux well, except for ID445. 

At 24 /»m, we find that the observed fluxes are typically 
higher than the model predictions, while the model flux of 
ID912 is consistent with the observed 24 /jm flux. We regard 
the model flux at 24 /im as the lower limit for the following 
reasons: 1) lack of AGN contribution in the model; and 2) 
the non-MW dust models are preferentially selected from the 
featureless SED at rest-frame UV, which may have a lower 
MIR emissivity than actual dust grains in submm galaxies. 
Specifically, the lack of AGN could be important for submm 
galaxies, in which the AGN activity to some extent is already 
known (Alexander et al. 2005). Recently, Daddi et al. (2007) 
found a galaxy population at z ~ 2 which shows a distinct 
excess of flux at 24 ^m, compared to that expected from the 
SFR estimated at other wavelengths. This MIR excess is 
attributed to Compton-thick AGNs, as a result of stacking 
analysis of deep X-ray images. Submm-bright NIRGs show- 
ing a clear excess of flux at 24 /im (ID300 and 1390), might 
be such MIR-excess galaxies. 

Submm-bright NIRGs are found to be massive, with the 
stellar masses of 5 x 10 10 -10 11 M Q . The derived bolomet- 
ric luminosity of 3 x 10 12 -10 13 L Q (excluding ID445) corre- 
sponds to the star formation rate of ~ 300-1000 M0yr _1 
for the adopted top-heavy IMF. For the Salpeter IMF, the 
stellar mass and the SFR could be even a few times higher 
than we derived. 

Following Takagi et al. (2004), we predict the present- 
day colours and absolute magnitudes of submm-bright 
NIRGs. Since we assume that the effects of star formation 
after the observed epoch is negligible, the derived present- 
day colours and luminosities are both lower limits. These 
lower limits could be close to reasonable values for ID300 
and 1390, which seem to be well evolved starbursts with 
t/to > 2 and M sta r > Af gas . For these galaxies, we predict ab- 
solute V-band magnitudes of My = —20.7 to —20.1, while 
submm-faint BzKs would have My = —19.4 on average. 
Since we use a top-heavy IMF, this magnitude is ~ 1 mag 
fainter than those by the Salpeter IMF. We find that the 
predicted rest-frame U — V and My are consistent with the 
colour-magnitude relation of elliptical galaxies. This is not 
the case for the Salpeter IMF (see also Takagi et al. 2004). 



6.2.3 Colour evolution of the SBURT model 

In Section 4, we mentioned the possibility that submm 
galaxies have a much larger overlap with BzA'-selected 
galaxies, compared with EROs and DRGs. In Figure 10, 
we show the BzK, R — K, and J — K colours as a function 
of starburst age at z = 2, i.e. a typical redshift for submm 
galaxies. We choose O = 0.7-1.4, which results in a good 
match with the observed SED variation of submm galax- 
ies for the SMC dust model (see Takagi & Pearson 2005). 
The SBURT models satisfy BzK > - 0.2 for a wide range 
of model parameters, i.e. starburst age and compactness of 
the starburst region (or optical depth). On the other hand, 
R—K > 3.35 is only satisfied with old models having t/to ^ 3 
for a wide range of O. This is also true for DRGs, except 
for a small fraction of models. Therefore, in the SBURT 
model, SEDs of stellar populations need to be intrinsically 
red to reproduce the colours of EROs and DRGs. Thus, 
this model predicts the largest overlap between BzKs and 
submm galaxies, among NIRGs. If we assume the limit on 
starburst age of t/to ~ 6 (Takagi et al. 2004), we expect the 
number of BzKs in submm galaxies to be about twice that 
of EROs, given that the selection of BzKs and EROs are 
not sensitive to the compactness O. This prediction could be 
confirmed by using larger samples of submm-bright NIRGs. 

6.2.4 Why are some BzKs bright in submm? 

In this study, we found that only a handful of NIRGs are 
bright in the submm. So, are submm-bright ones experienc- 
ing a special luminous evolutionary phase or simply more 
massive than submm-faint ones? If the former is the case, 
the duty cycle of submm-bright BzKs may be estimated 
from the fraction of submm-bright NIRGs. In the following 
discussion, we focus only on BzKs, since EROs and DRGs 
could be contaminated by passively evolving galaxies. We 
found that only 3 % of BzKs are bright in submm. This in- 
dicates that the duty cycle of submm-bright BzKs is only 
<10 % of the star- forming phase, even if we consider the in- 
completeness of the SHADES survey (Coppin et al. 2006). 
Thus, we might be observing a very sharp peak in the star 
formation activity of BzKs (e.g. see also Dannerbauer et al. 
2006). If this is the case, the SEDs of submm-bright BzKs 
may be systematically different from those of submm-faint 
ones. 

In order to address this question, we derived the aver- 
age SED of submm-faint BzKs with 1.5 < (B — z) < 2.5, i.e. 
having similar colour to submm-bright ones. In Figure 11, we 
show the average SED of 65 submm-faint BzKs with 1.5 < 
(B — z) < 2.5. At 850 /im, we adopted the average flux of 
BzKs derived from our stacking analysis. The average flux 
of submm-faint BzKs in the radio were also derived with a 
stacking analysis, and found to be /1.4GHz = 5.7± 1.0 fJ,Jy. A 
representative model 6 for the average SED was sought with 
the same method used for individual submm-bright galax- 
ies. For the flux error at optical-NIR bands, we adopted the 
standard deviation of the sample at each photometric band. 
Since this error is rather large, we fixed the redshift to the 

6 We call this model not the 'best' but 'representative', since the 
variance of the average SED is too large to specify one particular 
model as the best one. 
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average redshift of BzKs, (z) =1.9 (Daddi et al. 2004), in 
order to constrain the parameter space. We show the model 
thus obtained in Figure 11. The model parameters are tab- 
ulated in Table 5. 

Although the uncertainty is large, the average SED thus 
derived is not very much different from those of submm- 
bright BzKs, and reasonably explained by the the SBURT 
model, i.e. a starburst model, not mild evolutionary mod- 
els for quiescent galaxies. Therefore, from the SED analysis, 
we found no clear signatures that the evolutionary phase 
of submm-bright BzKs is substantially different from that 
of submm-faint ones. Compared to submm-bright galaxies, 
a difference in the model parameters may be found in the 
mass scale, i.e. submm-faint ones are less massive. Since the 
representative model is an old model of a starburst, the es- 
timated mass is close to the upper limit, owing to a higher 
mass-to- light ratio compared with younger models. From an 
oldest model which reproduces the average SED, we derive 
the upper limit on the stellar mass as M star <5 x 10 10 M . 
Compared to the oldest submm-bright NIRG in our sample, 
ID1390 (excluding ID445), the stellar mass of submm-faint 
BzKs could be less than half of ID1390. This may suggest 
that submm-bright and faint BzKs evolve into galaxies with 
different stellar masses, rather than being galaxies of a simi- 
lar stellar mass but in different evolutionary stages. We how- 
ever caution that the sample size of submm-bright BzKs are 
too small to derive firm conclusions on their average physical 
properties, such as stellar mass. We need a larger sample of 
submm-bright BzKs, and spectroscopic redshifts for more 
secure analyses. Also note that the mass estimate of the 
SBURT model may suffer from systematic effects, owing to 
the assumed simple star/dust geometry. 



7 SUMMARY 

We have investigated the submm properties of the following 
classes of NIR-selected galaxies: BzK-selected star-forming 
galaxies; DRGs; and EROs. We utilised a 93 arcmin 2 sub- 
region of the SHADES SXDF 850fim map which with has 
already been imaged in the NIR with the SIRIUS camera 
on the UH2.2m telescope. 

Using two SCUBA 850-^m maps (the SHADES B-map 
and D-map) produced by two different groups within the 
SHADES consortium, we detected 6 NIRGs above 3a. Four 
submm-detected NIRGs out of six are also detected both at 
24 /im and in the radio. This suggests that these 4 submm- 
detected NIRGs are genuine submm-bright galaxies. These 
submm-bright NIRGs are all BzK-selected galaxies, except 
for ID1390 whose BzK colour is however close to the se- 
lection boundary. In other words, no EROs and DRGs are 
found to be submm-bright if they are clearly non-BzKs. 
We made a rough estimate of the number of BzKs in radio- 
detected submm galaxies, assuming that submm galaxies at 
1.4 < z < 2.5 satisfy the selection criteria of BzKs. Al- 
though the sample is small, this estimate is consistent with 
our result. This may indicate that most submm galaxies at 
1.4 <z < 2.5 could be BzKs. 

Two submm-detected NIRGs satisfy all the selection 
criteria we adopted, i.e. they are extremely red BzK- 
selected galaxies. Although these extremely red BzKs are 
rare, the fraction of submm galaxies in them could be 



high (up to 20%), compared with the other colour-selected 
optical-NIR galaxy populations. 

We performed stacking analyses with our SCUBA 850- 
/im maps, in order to derive the average flux of submm-faint 
NIRGs. We derived 0.52±0.19 (0.64±0.16) mjy and 0.53± 
0.16 (0.50 ± 0.13) mjy from the B-map (D-map) for BzKs 
and EROs, respectively, while we found no significant signal 
from DRGs in either map. The contribution from BzKs and 
EROs to the EBL at 850 /jm is about 10-15%. Focusing on 
the EBL only from BzKs, we found that > 30 % of the EBL 
from BzKs is resolved in our SCUBA map. This is higher 
than that for EROs and submm sources as a whole. This 
might be expected if the majority of the EBL originates 
from submm-undetected galaxies at z < 1.5, as suggested 
by Wang et al. (2006). For galaxies at z ~ 2, the fraction of 
submm flux in resolved sources could be higher than that in 
low- z galaxies. 

We have also fitted SED models of starbursts to each 
of the submm-bright NIRG (mostly BzKs) and to an aver- 
age SED of submm-faint BzKs derived from galaxies with 
1.5 < (B — z) < 2.5, i.e. similar colours to the submm- 
bright ones. Submm-bright NIRGs are found to lie at the 
typical redshifts of submm galaxies and have stellar masses 
of 5 x 10 10 -10 n M Q with a Salpeter-like, but slightly top- 
heavy IMF and bolometric luminosity of 3 x 10 12 -10 13 L Q . 
From the average SED of submm-faint BzKs, we found no 
clear signature that the evolutionary phase of submm-bright 
BzKs are substantially different from that of submm-faint 
ones, as suggested by the small number of submm-bright 
BzKs. On the other hand, submm-faint BzKs are likely to 
be less massive, with the stellar mass below - 5 x 10 10 M . 

The results presented here are clearly limited by small 
sample size. Nevertheless, our study can still be consid- 
ered useful for investigating the physical relationships be- 
tween NIR-selected and submm-selected massive galaxies. A 
large sample of submm-bright NIRGs from currently ongo- 
ing and future muti-wavelength surveys, including SHADES 
and UKIDSS will play an important role on the study of 
massive galaxy formation and evolution. 
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APPENDIX A: ESTIMATE ON THE NUMBER 
OF CHANCE DETECTIONS IN OUR SCUBA 
MAPS 

The large beam size of the JCMT (14.7" FHWM at 850 fj,m) 
could lead to some chance detections of spurious >3<r peaks 
at the position of NIR-selected galaxies. We statistically 
evaluate the expected number of >3<r chance detections for 
a given class of galaxies as follows. The probability of ran- 
domly finding an object at >3<r in the 850-/im map may be 
given by the ratio of the number of pixels with >3<r to the 
total number of pixels in the region under consideration. For 
the SHADES/SIRIUS A^-band area, we found this proba- 
bility as p(>3a) = 0.45 x 10" 2 for the SHADES B-map and 
0.85 x 10~ 2 for the D -map . If we randomly distribute n ob- 
jects in the 850-/xm map, the number of objects spuriously 
detected with >3 a would follow a Poisson distribution with 
the parameter of /i = np(>3u). For 307 NIRGs, we derive 
\x = 1.4 and 2.6 for the B- and D-map respectively, giving 
the mean number of >3 a detections by chance. This means 
that we must confirm whether submm-detected NIRGs are 

7 The difference in p(>3cr) between the maps may be explained 
by the difference in the adopted pixel size. 
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Table 1. Statistics of NIR-selected galaxies (NIRGs) 



Population 


Area" 


Number 


Surface density 




K B or J 




[arcmin - 2 ] 


-ff a -detected object 


K s 


1308 (992) b 


11.5 


Star 


K s 


95 (64) 


0.83 


NIRGs 


J 


307 (245) 


3.14 


ERO 


Kg 


249 (201) 


2.18 


DRG 


J 


84 (67) 


1.09 


ERO and non-DRG 


J 


105 (84) 


1.36 


DRG and non-ERO 


J 


38 (28) 


0.49 


ERO and DRG 


J 


46 (39) 


0.60 


BzK 


K s 


168 (132) 


1.47 


BzK and ERO and non-DRG 


J 


17 (14) 


0.22 


BzK and DRG and non-ERO 


J 


14 (8) 


0.18 


BzK and ERO and DRG 


J 


19 (16) 


0.25 



a) Image (K B - or J-band area) for which the number of objects are derived. 



b) Total number of objects, with the number within the SHADES area given 
in parentheses. 



Table 2. Submm properties of submm- bright NIR-selected galaxies 



Source 


Map a 


S/N b 


nc 

°850Mm 


cd 

°850Mm 


■S'24fim 


<Sl.4GHz 


SHADES ID e 








[arcsec] 


[mjy] 


[MJy] 


[/uJy] 


[SXDF] 


300 


B,D 


3.1 


4.5 


5.7± 2.1 


546±27 


27.8± 7.0 


850.30 


445 


B 


3.3 


2.9 


5.6± 2.1 


357±24 


250.± 7.4 


850.27 


912 


B,D 


1.1 


2.1 


4.4± 1.8 


511±27 


145. ± 7.7 


850.4 


1390 


B,D 


3.4 


1.8 


4.0± 2.8 


446±29 


29.6± 7.5 


(850.62)-^ 


Unconfirmed detections 


718 


B,D 


3.3 


2.3 


4.0± 2.1 


< 450 


< 35 


850.70 


1133 


D 


3.0 


6.6 


3.0± 2.1 


< 450 


< 40 


850.77 



Notes: The upper limits correspond to 5 cr. a) SHADES map in which NIRGs are detected 

b) S/N at the position of NIRGs in the SHADES map. We adopt values from the SHADES 
B-map if NIRGs are detected in both B and D maps, c) The angular separation between 
submm and NIR position. Submm positions are taken from the SHADES source catalogue. 
For ID1390, we adopt the submm position from the SHADES B-map. d) Deboosted 850-/im 
flux densities from the SHADES source catalogue, using the algorithm of Coppin et al. 
(2005). e) Source names in the SHADES source catalogue, f) Source name from a preliminary 
SHADES source catalogue. 



genuine submm emitters or not by other means. Here we re- 
quire a detection at 24 /im or radio to identify submm-bright 
NIRGs, along with the detection in the SCUBA map. For 
example, only 15 NIRGs are found to be detected at 24 /im. 
This results in fj, = 0.067 and 0.13 for 24 /mi-detected NIRGs 
in the B-map and D-map, respectively. Thus, we expect no 
chance detections in the SCUBA map if NIRGs are detected 
at 24 fim. 
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Table 3. Optical - NIR properties of submm-bright NIR-selected galaxies 



Source 


R.A. 


Dec. 


K B 


R-K B 


J-K B 


B-z' 


z' -K s 




Notes 








[J2000] 


[J2000] 






[AB mag] 






BzK 


ERO 


DRG 


24 iim 


Radio 


300 


2 17 40.01 


-05 01 15.7 


21.14 


3.43 


1.56 


2.05 


2.61 


Y 


Y 


Y 


Y 


Y 


445 


2 18 07.92 


-05 01 45.7 


22.00 


4.39 


> 1.40 


1.75 


3.21 


Y 


Y 


Y 


Y 


Y 


912 


2 17 38.67 


-05 03 39.4 


21.64 


2.59 




1.20 


2.05 


Y 


N 




Y 


Y 


1390 


2 18 07.74 


-05 06 10.5 


21.13 


3.67 


2.02 


2.52 


2.56 


(N) a 


Y 


Y 


Y 


Y 



Notes: A 'Y' ('N') in the last five columns indicate that the galaxy is (not) Bzff-selected galaxy, ERO, DRG, 24 ^tm-detected 



and radio-detected, a) After the colour correction on B — z and z — K, this source becomes a non-Bz-ff-selected galaxy, as 
shown in Figure 7. This correction is adopted to match the observed colour sequence of stars to that of Daddi et al. (2004). 



Table 4. Summary of stacking analysis at 850 ^tm 



Class 


Number" 


(Ss50lJ.ni) b 


K-S Prob. c 


% Resolved flux d 


E.B.L. e 






[mjy] 


[%] 


[%] 


[Jydeg- 2 ] 


BzKs 


112 


0.52 ± 0.19 


5.91 


32 ± 18 


3.8 ± 1.2 


EROs 


178 


0.53 ± 0.16 


1.70 


18 ±9 


5.1 ± 1.5 


DRGs 


56 


0.30 ±0.28 


45.91 







a) The number of objects used for the stacking analysis, b) The average 850-/^m 



flux of non-detected objects, which arc > 7" away from individual sources, c) The 
probability from the Kolmogorov-Smirnov test, that the flux distribution of objects 
are derived from the same sample, d) The fraction of the flux density in individually 
detected sources, e) The extragalactic background light at 850 (im from each class of 
objects. 



Table 5. Summary of the SED fitting with the SBURT model 



Source 


X 2 /v a 


^•phot 


Age 6 


e c 


Ext. d 


log M s tar 


log M gas 




A v 


log SFR 


All 


(U-V)f 








[Gyr] 






[M ] 


[M Q ] 


[L©1 


[mag] 


[M Q yr- 1 ] 


[mag] 


[mag] 


300 


0.91 


2.8 


0.2 


2.0 


SMC 


10.8 


10.7 


12.7 


1.3 


2.7 


-20.1 


1.3 


445 


1.47 


2.1 


0.5 


1.0 


LMC 


10.9 


10.0 


12.1 


2.3 


2.0 


-20.3 


1.5 


912 


0.47 


2.1 


0.07 


1.4 


LMC 


10.7 


11.1 


13.1 


3.1 


3.2 


-20.1 


0.9 


1390 


2.45 


2.1 


0.4 


1.4 


LMC 


11.1 


10.4 


12.5 


1.5 


2.1 


-20.7 


1.4 


(BzKs) 3 




1.9 h 


0.4 


1.6 


LMC 


10.6 


9.8 


11.9 


1.1 


1.8 


-19.4 


1.4 



a) x 2 divided by the degree of freedom v. b) Starburst age with the evolutionary time scale of to = 0.1 Gyr. 

c) Compactness parameter of starburst region, d) Extinction curve used, e) Bolometric luminosity, f) Predicted 
present-day K-band magnitude and U — V (Vega), assuming the passive evolution after the observed epoch, g) 
Model parameters for an average SED of submm-faint BzKs with 1.5 < (B — z) < 2.5 h) Assumed redshift for the 
SED fitting. 
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Figure 1. J — K B versus R — K B for A" s -band detected objects in the SIRIUS J-band area. The selection criteria for EROs and DRGs 
are indicated with a vertical line and a horizontal line, respectively. 
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K [AB] 

Figure 2. Bzi^" colour vs. _R" S magnitude for _ff s -band detected objects. Here the z' — K B and B — z' colours have not been colour- 
corrected in order to show the original photometric errors. The BzK selection criterion proposed by Daddi ct al. (2004) corresponds to 
[z' — K B ) — (B — z') = 0.2 for the adopted colour correction (see text). 




ARC SECONDS 10 5 -5 -10 

CENTER: R.A. 02 18 7.94 DEC -05 01 45.6 J2000 arc seconds 

Figure 3. [Left column] Contours of SCUBA 850-^tm for submm-bright NIR-selected galaxies. The corresponding S/N ratios at 850 fim 
are indicated at each contour line. The underlying images are fC s -band with 20" on a side. The names of the sources are indicated at the 
top of each image. [Right column] Contour maps of radio (1.4 GHz) on negative images at ii-band for the same regions as shown in the 
left panels. Contour lines are shown at —2, 2, 3, 4, ... ,10, 20, ... 100 Xcr. Positive (negative) contours are indicated with solid (dashed) 
lines. 
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CENTER: R.A. 02 17 38.67 DEC -05 03 39.4 J2000 
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CENTER: R.A. 02 18 7.77 DEC -05 06 10.7 J2000 

Figure 3. - continued. 
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Figure 4. Histograms of the 850- /xm flux at the position of BzKs (shaded; right axis), compared to the histograms for the unmasked 
regions of the SCUBA map as a whole (open; left axis). 
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Figure 5. Same as Figure 4, but for EROs. 
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Figure 6. Same as Figure 4, but for DRGs. 




Figure 7. The BzK diagram for if s -detected sources in the SXDF/SIRIUS field. Solid circles indicate submm-bright NIRGs. The 
photometric errors include those of colour corrections on B — z and z — K. Both EROs and DRGs are indicated with small circles. Large 
thick/thin crosses are for NIR-selected/i^s-detected galaxies which are detected at 24 fim with SWIRE. Small pluses indicate stars from 
Daddi's catalogue. We depict the B,zii"-selection criterion with a solid line. Dashed and dotted lines are boundaries for selecting passively 
evolving galaxies and stars, respectively. 
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Figure 8. Solid and dotted linos indicate the best-fitting SBURT model without and with an extended radio component of the SED, 
respectively. The data points used for the SED fitting are shown as solid circles. Arrows indicate the 5<r upper limits. See Section 6.2.1 
for the definition of the fitting parameters. 
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Figure 9. Contour maps of Ax 2 from the SED fitting in a plane of starburst age (with to = 0.1 Gyr) and redshift. The contours are 
depicted at A\ 2 = 2.71, 6.63, having the probabilities of 90 and 99 % when projected on to each axis. We note that there are no local 
minima at z < 1 and z > 4. 
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Starburst age (t/t ) 

Figure 10. BzK, R — K and J — K colours of the SBURT model at z = 2 as a function of starburst age, which show the range of model 
parameters to satisfy BzK, DRG, and ERO selection criteria. Model parameters for a given line symbol are shown in the legend box. 
Horizontal lines indicate the colour boundary for BzKs, DRGs and EROs. 
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Figure 11. Average SED of BzKs with 1.5 < (B — z') < 2.5 and a representative SBURT model for the average SED. Symbols are the 
same as in Figure 8. 



